Probing entanglement in a 2D hard-core Bose–Hubbard lattice

Entanglement and its propagation are central to understanding many physical properties of quantum systems1–3. Notably, within closed quantum many-body systems, entanglement is believed to yield emergent thermodynamic behaviour4–7. However, a universal understanding remains challenging owing to the non-integrability and computational intractability of most large-scale quantum systems. Quantum hardware platforms provide a means to study the formation and scaling of entanglement in interacting many-body systems8–14. Here we use a controllable 4 × 4 array of superconducting qubits to emulate a 2D hard-core Bose–Hubbard (HCBH) lattice. We generate superposition states by simultaneously driving all lattice sites and extract correlation lengths and entanglement entropy across its many-body energy spectrum. We observe volume-law entanglement scaling for states at the centre of the spectrum and a crossover to the onset of area-law scaling near its edges.


I. INTRODUCTION
Entanglement is a uniquely quantum property that underpins descriptions of interacting quantum systems as statistical ensembles [1][2][3][4].Within closed many-body quantum systems, entanglement amongst constituent subsystems introduces uncertainty to their individual states, even when the full system is in a pure state [5,6].For this reason, entropy measures are commonly used to quantify quantum entanglement in interacting manybody systems and have been directly probed in a number of platforms [7][8][9][10][11][12][13]. The study of entanglement in interacting many-body quantum systems is central to the understanding of a range of physical phenomena in condensed-matter systems [5], quantum gravity [14,15], and quantum circuits [16].The scaling of the entanglement entropy with subsystem size provides insight into classifying phases of quantum matter [17][18][19] and the feasibility of numerically simulating their time dynamics [6].
In a closed system, the bipartite entanglement entropy of a subsystem quantifies the amount of entanglement between the subsystem and the remainder of the system.For certain many-body states, such as the ground state of 1-dimensional local Hamiltonians [20], the entanglement entropy is proportional to the size of the boundary between a subsystem and the remaining system; this boundary is referred to as the area of a subsystem.Such states are said to have area-law entanglement scaling [6].
For other states, entanglement entropy increases proportionally to the bulk size (volume) of a subsystem, a behavior referred to as volume-law scaling.To characterize the entanglement scaling in a lattice of quantum objects, such as natural or artificial atoms, we consider the area and volume entanglement entropy per lattice site, represented by s A and s V , respectively.Disregarding logarithmic corrections, which are theoretically expected in certian contexts [21,22], for a given state, the entanglement entropy S(ρ X ) of a subsystem X can be expressed using the ansatz where A X is the subsystem area and V X is the subsystem volume (see Fig. 1a for an example).The ratio of volume to area entropy per site, s V /s A , quantifies the extent to which that state obeys area-or volume-law entanglement scaling [23].Quantum states with area-law entanglement scaling have local correlations, whereas systems obeying volume-law entanglement scaling contain correlations that extend throughout the system and are therefore more challenging to study using classical numerical methods [24,25].
In this work, we study the entanglement scaling of states residing in different energetic regions of the twodimensional (2D) hard-core Bose-Hubbard (HCBH) lattice.The Bose-Hubbard model is particle-number conserving, allowing its energy spectrum to be partitioned into sectors with definite particle number n.The "hard-core" condition arises from strong on-site particleparticle interactions that mandate that each lattice site may be occupied by at most a single particle (Fig. 1b).In 2D, the HCBH model is non-integrable and may ex- hibit eigenstate thermalization [2,26,27] and quantum information scrambling [28].Highly excited many-body states-states residing near the center of the energy spectrum (energy E ≈ 0; orange oval in Fig. 1c)-are expected to exhibit volume-law scaling of the entanglement entropy, following the Page curve [14] (orange line in Fig. 1d).For subsystems smaller than half the size of the system, the entropy grows linearly in subsystem size.
The entropy is maximized when the subsystem comprises half of the entire system and decreases as the subsystem size is further increased.In contrast, the entanglement entropy of states residing near the edges of the HCBH energy spectrum does not follow the Page curve, but rather, grows less rapidly with subsystem volume (exemplified by the teal line in Fig. 1d showing the entropy of the eigenstate highlighted by a teal oval in Fig. 1c).For states with intermediate energies, a crossover is expected from area-law scaling at the edges of the energy spectrum to volume-law scaling at its center (Fig. 1e) [22,23].Although Fig. 1e illustrates the crossover for the n = 8 particle-number manifold, the crossover similarly occurs in other manifolds (see numerical simulations in Fig. S12 of the Supplementary Material).
The crossover from area-to volume-law entanglement scaling can be observed by studying the exact eigenstates of the HCBH model.However, preparing a specific eigenstate generally requires a deep quantum circuit [29] or a slow adiabatic evolution [30].Alternatively, we can explore the behavior of the entanglement entropy by preparing superpositions of eigenstates across multiple particle-number manifolds [23].Here, we prepare such superposition states by simultaneously and weakly driving 16 superconducting qubits arranged in a 4×4 lattice.By varying the detuning of the drive frequency from the lattice frequency, we generate superposition states occupying different regions of the HCBH energy spectrum.Measurements of correlation lengths and entanglement entropies indicate volume-law entanglement scaling for states prepared at the center of the spectrum and a crossover to the onset of area-law entanglement scaling for states prepared at the edges.
In this experiment, we use a 2D lattice of 16 capacitively-coupled superconducting transmon qubits [37] fabricated using a flip-chip process [38], as depicted in Fig. 1f.The qubits are located on a qubit tier (Fig. 1g), and the readout and control lines are located on a separate interposer tier (Fig. 1h, see Supplementary Material for device details).The circuit emulates the 2D Bose-Hubbard model, described by the Hamiltonian where â † i (â i ) is the creation (annihilation) operator for qubit excitations at site i, ni = â † i âi is the respective excitation number operator, and qubit excitations correspond to particles in the Bose-Hubbard lattice.The first term represents the site energies ϵ i = ω i − ω r , which are given by the transmon transition frequencies ω i with respect to a frame rotating at frequency ω r .The second term accounts for the on-site interaction with strength U i arising from the anharmonicity of transmon qubit i, representing the energy cost for two particles to occupy the same site with an average strength U/2π = −218 (6) MHz.The final term of the Hamiltonian describes the particle exchange interaction (within a rotating wave approximation) between neighboring lattice sites with strength J ij , realized by capacitively coupling adjacent qubits with an average strength of J/2π = 5.9(4) MHz at qubit frequency ω/2π = 4.5 GHz.While the coupling strengths between qubits are fixed, we switch particle exchange off for state preparation and readout by detuning the qubits to different frequencies (inset in Fig. 2a).Our system features site-resolved, multiplexed single-shot dispersive qubit readout [39] with an average qubit state assignment fidelity of 93%, which, together with site-selective control pulses, allows us to perform simultaneous tomographic measurements of the qubit states.
In our system, on-site interactions are much stronger than exchange interactions, J ≪ |U |.By restricting each site to two levels, and mapping the bosonic operators to qubit operators, we transform the Hamiltonian in Eq. 2 to the hard-core limit: where σ+ i (σ − i ) is the raising (lowering) operator for a qubit at site i, and σz i is the Pauli-Z operator.The energy relaxation rate Γ 1 and dephasing rate Γ ϕ in our lattice are small compared to the particle exchange rate, with Γ 1 ≈ 10 −3 J and Γ ϕ ≈ 10 −2 J, which enables us to prepare quantum many-body states and probe their time dynamics faithfully.

III. COHERENT-LIKE STATES
We generate a superposition of many eigenstates of the lattice, which we refer to as a coherent-like state, by simultaneously driving the qubits via a common control line.Applied to the lattice initialized with no excitations, the drive acts as a displacement operation of the Fock basis defined by the number of excitations in the lattice [23], hence, it will result in a state analogous to coherent states of light.The common drive is applied to the system via the feedlines coupled to each qubit through their respective readout resonators.Selecting the rotating frame frequency ω r to be the frequency of the drive, the Hamiltonian of the driven lattice is where δ = ω r − ω com is the detuning between the drive and the qubit frequencies (all qubits are biased onresonance at ω com ).The drive strength Ω can be tuned by varying the amplitude of the applied drive pulse.The common drive couples independently to each qubit with a complex coefficient α j that depends on the geometric circuit parameters of the lattice.
We first study the time dynamics of the average number of excitations ⟨n⟩ in the lattice under a resonant drive (δ = 0) in Fig. 2a.The driven lattice reaches a steady state of half-filling, with an average particle number ⟨n⟩ = 8, after driving the lattice for time t ≈ 10/J.Once in steady-state, the drive adds and removes excitations coherently from the system at the same rate.Hamiltonian parameters J ij and α j were characterized through a procedure described in Section S4 of the Supplementary Material, and the excellent agreement of numerical simulations of time evolution under Eq. 4 with experimental data confirms their accuracy.In Fig. 2b, we report the discrete probability distribution of measuring a different number of excitations in the lattice at three different times.The probability of a particular excitation number approximately follows a Poisson distribution for a weak drive Ω = J/2 (blue stars with dashed lines as guides to the eye), signifying that the quantum state is in a coherent-like superposition of excitation-number states.
The coherent-like state comprises a swath of the HCBH energy spectrum that depends on the drive detuning.Driving the lattice on resonance prepares a coherentlike state at the center of the HCBH energy spectrum (Fig. 2c).By varying the drive detuning, we can generate states that are a superposition of the eigenstates closer to the edge of the energy band (Fig. 2d,2e).The standard deviation in the energy of the coherent-like state depends on the strength of the drive: a stronger drive increases the bandwidth of populated energy eigenstates within each particle number subspace.Therefore, to probe the entanglement properties across the lattice energy spectrum, we choose a relatively weak drive with strength Ω = J/2 for the remainder of the main text.We choose a drive dura- tion t = 10/J, which is short compared to the timescale of decoherence yet long enough to allow the coherent-like state to reach its steady-state distribution (which occurs at roughly t = 6/J according to simulations shown in Section S15 of the Supplementary Material).

IV. CORRELATION LENGTHS
A system of interacting particles exhibits correlations between its constituent subsystems.The hard-core Bose-Hubbard Hamiltonian is equivalent to an XY Hamiltonian, where quantum order is reflected by the transverse correlations [40].In order to quantify transverse correlations for coherent-like states generated with detuning δ, we make single-shot measurements of identically prepared systems and then calculate the 2-point correlators C x i,j ≡ ⟨σ x i σx j ⟩ − ⟨σ x i ⟩⟨σ x j ⟩ between different qubit pairs.In Fig. 2f, we show the magnitude-squared 2point correlator values |C x i,j | 2 , averaged over qubits at the same Manhattan (1-norm) distance M .When δ/J is small, generating a superposition state occupying the center of the energy band, |C x i,j | 2 becomes small for all M , matching the expectation for states with volume-law entanglement scaling.The 2-point correlator has an upper-bound set by the mutual information between the sites I(i : j) = S(ρ i ) + S(ρ j ) − S(ρ ij ).In a volume-law state, the entropy of small subsystems is equal to their volume [14], hence we expect the mutual information and, in turn, the correlator between any qubit pair will vanish.
Using the 2-point correlators, we extract the correlation length ξ x by fitting |C x i,j | 2 ∝ exp(−M/ξ x ) (Fig. 2g).The correlation length quantifies the dependence of correlations on the distance between the subsystems within our system.As we increase the magnitude of the drive detuning δ, skewing the superposition state to the edge of the energy spectrum, we generally find that the correlation length increases.The states prepared with −J/2 < ∼ δ < ∼ J/2, however, follow the opposite trend, and the extracted correlation length diverges around δ = 0.In this regime where the coherent-like state occupies the center of the energy spectrum, the 2-point correlators asymptotically reach zero, so the extracted 2-point correlation length loses meaning.We attribute the slight asymmetry of the observed correlation lengths about δ = 0 to a slight offset of the energy spectrum of our lattice towards positive energies, which, as shown in Section S5 of the Supplementary Material, is a consequence of nextnearest neighbor exchange interactions.

V. ENTANGLEMENT SCALING BEHAVIOR
In order to study the entanglement scaling of the coherent-like states, we reconstruct the density matrices of 163 unique subsystems, listed in Section S9 of the Supplementary Material, via a complete set of tomography measurements.The measured subsystems contain up to 6 qubits (see Fig. 3a for examples).We quantify the entanglement entropy between subsystem X and the remaining system through the second Rényi entropy [41] where ρ X is the reduced density matrix of a subsystem X of the entire quantum system ρ.
We first study the scaling of the entanglement entropy with the subsystem volume V X , defined as the number of qubits within the subsystem, for coherent-like states prepared with different drive detunings δ (Fig. 3b).For coherent-like states prepared with δ = 0, we observe nearly maximal scaling of the entropy with subsystem size S 2 (ρ X ) ≈ V X , whereas with increasing δ the entropy grows more slowly with subsystem size (Fig. 3c).
There is excellent agreement between the expected and the measured entropy for subsystems of volume ≤ 4, yet there is a discrepancy for the largest subsystems where the coherent-like states are prepared near the center of the spectrum (see Fig. 3c).The discrepancy arises, in large part, from having a finite number of measurement samples when reconstructing subsystem density matrices.Tomographic reconstruction of the state of a subsystem requires determining the distribution of measurement outcomes for each Pauli string.When the system has volume-law entanglement, the measurement outcome distribution becomes more uniform.Density matrix reconstruction therefore becomes more sensitive to sampling error.This becomes especially relevant as the subsystem size increases.Therefore, for the largest subsystems, the entropy of the reconstructed density matrix will be smaller than the actual entropy of the subsystem, as seen for V = 5, 6 in Fig. 3c.We simultaneously measured all subsystems by taking 2,000 measurement samples for 3 6 Pauli strings, as visualized in Fig. S15 of the Supplementary Material.For subsystems of volume V , we can therefore extract 2,000×3 6−V measurement samples-sufficient for V = 1-4, but less so for V = 5, 6.In hindsight, after the experiment had concluded, we realized that we could have obtained better agreement for V = 5, 6 if we had used 20,000 samples (open diamonds in Fig. 3c), which would have been straightforward to implement experimentally.In Section S11 of the Supplementary Material, we discuss this sampling dependence in more detail and show, using Monte Carlo simulations, that a more accurate reconstruction of highly entangled states follows from larger numbers of measurement samples of each Pauli string.
We next determine the scaling of entanglement entropy with subsystem volume s V and area s A per Eq. 1.Using the Rényi entropies of the density matrices reconstructed from experimental data, we extract s V and s A by measuring the rate of change of S 2 (ρ X ) with V X and the subsystem area A X , respectively, with the other held fixed.Here, A X is defined as the number of nearestneighbor bonds intersecting the boundary of the subsystem X.The linear fitting procedure used to determine s V,A is detailed in Section S10 of the Supplementary Material.In Fig. 3d we observe that as the magnitude of δ becomes larger, s V decreases and s A increases.While extraction of s V is reliable at all drive detunings, at small drive detuning values (−J < δ < J) the entanglement entropy does not exhibit a notable dependence on the area of the subsystem in our finite lattice, hence we are not able to reliably fit s A .By considering the geometric entropy ratio s V /s A we observe a change in the behavior of entanglement entropy within our system (Fig. 3e).The states prepared at the center of the energy spectrum with δ = 0 exhibit a strong volume-law scaling, and the states prepared closer to the edge of the energy spectrum show weaker volume-law scaling and increasing area-law scaling.
The entanglement spectrum, obtained via Schmidt decomposition, further quantifies the structure of entanglement across a bipartition of a closed quantum system [17].The quantum state |ψ⟩ can be represented as a sum of product states of the orthonormal Schmidt bases |k X ⟩ for a given subsystem X and |k X ⟩ for the remaining lattice [42]: where positive scalars λ k are the Schmidt coefficients with k λ 2 k = 1.The Schmidt coefficients form the entanglement spectrum and provide a proxy for the degree of entanglement between the two subsystems.For a subsystem X maximally entangled with the remaining lattice, all the Schmidt coefficients will have an equal value We obtain the entanglement spectrum for a bipartition of our lattice by diagonalizing the measured partial density matrix of a subsystem.In Fig. 4a, we study the entanglement formed within states prepared across the energy spectrum.We report the first 16 Schmidt coefficients squared λ 2 k , in decreasing order, for a subsystem highlighted in maroon and the remaining lattice.We observe that for states obeying volume-law entanglement scaling at the center of the spectrum, the variation in coefficient magnitudes is small compared to states closer to the edge of the spectrum, in close agreement with numerical simulation.To quantify this difference, in Fig. 4b we show the ratio of the largest and the k th largest Schmidt coefficient, λ 2 1 /λ 2 k , for k = 5, 10, 14, of coherent-like states prepared with different drive detunings δ.We observe that a small number of Schmidt states contain nearly all the weight of the decomposition for the area-law-like states, whereas, for the volume-law states, the Schmidt coefficients are roughly equal.This variation signals a change in the extent of the entanglement distribution across the system.
In Fig. 4c, we report the number of coefficients required to approximate the state of the lattice bipartition with accuracy 1 − ϵ = 0.999 (see Section S13 in Supplementary Material) for subsystems with volume V = 3, 4, 5.We find that the states at the edge of the energy spectrum can be accurately represented with fewer coefficients, while the number of coefficients needed to approximate states at the center of the band approaches the dimension of the Hilbert space of the subsystem, 2 V .

VI. CONCLUSION
In this work, we study the entanglement scaling properties of the 2D hard-core Bose-Hubbard model, emulated using a 16-qubit superconducting quantum processor.By simultaneously driving all qubits, we generate coherent superposition states that preferentially incorporate eigenstates from regions of the many-body energy spectrum that we tune between the center and edges.We probe the transverse quantum correlation lengths and the entanglement scaling behavior of the superposition states.We observe a crossover from volume-law scaling of entanglement entropy near the center of the band, marked by vanishing two-point correlators, to the onset of area-law entropy scaling at the edges of the energy band, accompanied by finite-range correlations.
A coherent-like state comprises eigenstates across a swath of the HCBH energy spectrum.Decreasing the drive amplitude will decrease the width of the coherentlike state.Because excitations (particles) are added more slowly to the lattice, the coherent-like state will need to evolve for more time before reaching steady-state.Therefore, coherent-like states can be prepared with narrower spectral width as the coherence of superconducting processors improves, providing finer energy resolution of the area-to-volume-law crossover.
A central outstanding question is when and how statistical ensembles arise from unitary time evolution [1-3, 43, 44].While an analytical relation between entanglement and thermodynamic entropy exists for integrable systems [45], an analogous insight for non-integrable systems could lead to a deeper understanding of the emergence of thermodynamic behavior in quantum systems.In recent years, random circuit protocols have emerged as a potential tool for addressing this question in a digital context [46].The protocol introduced in the present work may be viewed as an analog (i.e.continuous time evolution) counterpart to random circuit experiments, which can be utilized to generate highly entangled states rapidly.While tomography itself has exponential time complexity, our techniques are extensible to larger system sizes.In Section S12 of the Supplementary Material, we present numerical simulations indicating that the area-to-volume-law transition should remain clear in larger systems without increasing the volume of subsystems used for tomography.Because the maximum subsystem volume is fixed and the tomographic readout of all subsystems is simultaneous, our protocol has constant runtime as a function of system size.Therefore, as quantum technology moves to larger systems, but before full error correction is available, our protocol will elucidate emergent thermalization in a regime where results are classically intractable.
Finally, understanding the structure of entanglement within a quantum system determines the effective degrees of freedom that need to be considered to simulate the quantum states.Area-law states can generally be numerically simulated efficiently using tensor network methods [6,24,25], whereas the computational complexity of classically simulating volume-law states scales exponentially with system size.It is the latter complexity that underpins the promise of quantum computational advantage.In this work, we have introduced an extensible, hardware-efficient analysis of the scaling structure of entanglement that could serve as the basis for an algorithmagnostic means to verify the computational complexity of programs executed on large-scale quantum devices.
Note: During the preparation of this manuscript, we became aware of related studies in a 1D trapped-ion simulator [47].We house our superconducting qubit array within a BlueFors XLD-600 dilution refrigerator.In our setup, the mixing chamber (MXC) base temperature reaches as low as 14mK, however, we use a heater with a PID controller to stabilize the base temperature at 20mK.This mitigates temperature fluctuations during experiments.
We use a 24-port microwave package [1] to host the sample, consisting of a copper casing, a multilayer interposer, a shielding cavity in the package center, and a set of microwave SMP connectors.The signal and ground connections between the package and the sample are provided using superconducting aluminum wirebonds.
Semi-rigid coaxial cables with 0.086" outer conductor diameter and SMA connectorization transmit qubit control and readout signals between the mixing chamber stage and room temperature electronics.Cables connecting stages at different temperatures have low thermal conductivity to minimize the heat flow from warmer stages.We use KEYCOM ULT-05 cables between room temperature and the 4 K stage.All input lines (readout input, qubit control, and TWPA pump lines) use stainless-steel (SS-SS) cables between the 4 K and the mixing chamber stages.Readout output lines use low-loss superconducting niobium-titanium (NbTi) cables between the 4 K and mixing chamber stages to maximize the signal-to-noise ratio of the readout signal.We use hand-formable, non-magnetic, copper coaxial cables (EZ Form Cable, EZ-Flex.86-CU) between the mixing chamber stage and the sample package.
We attenuate and filter the signal at different stages of the dilution refrigerator to reduce the noise incident on the qubits.For the qubit control lines, we use a 20 dB attenuator (XMA Corporation) at the 4 K, a 1 dB attenuator at the still, and a 1 GHz low-pass filter (Mini-Circuits VLFG-1000+) at the mixing chamber stage.The 1 dB attenuator at the still was selected to thermalize the center conductor of the coaxial line while maintaining the ability to apply DC current through the line without causing excessive heating.The VLFG-1000+ low-pass filter reduces the noise at higher frequencies while allowing an RF drive to reach the qubit.The filter response in the range between 4 − 5 GHz varies by only 3 dB, approximately providing 40 dB of attenuation.Our readout input lines contain a 20 dB attenuator at 4 K, a 10 dB attenuator at the still, and 40 dB of attenuation at the mixing chamber stage, followed by a 12 GHz low-pass filter(RLC F-30-12.4-2).This attenuation and filtering scheme mitigates the thermal-photon population in the resonators [2], which limits the qubit coherence times due to photon shot noise.
The readout circuit in the samples used in our experiment is set up for measuring in reflection mode.The resonator probe tone is coupled to the sample via a circulator (LNF-CIC4 12A).The reflected readout signal gets routed to the readout output chain using the third port of the circulator.The signal is amplified at the mixing chamber stage using a near-quantum-limited Josephson traveling-wave parametric amplifier (TWPA) [3] fabricated at MIT Lincoln Laboratory and mounted inside the µ-metal shield.The TWPA pump tone is combined with the readout signal using a directional coupler (Marki C20-0116).The signal then travels through two isolators (Quinstar 0XE89) followed by a 12 GHz low-pass filter (RLC F-30-12.4-2) and a 3 GHz high-pass filter (RLC F-19704) to prevent noise from coupling back to the sample through the readout chain.Afterward, the signal is once again amplified using a low-noise highelectron-mobility transistor (HEMT) amplifier at the 4 K stage before being transmitted outside the cryostat to room temperature measurement electronics.
We use a microwave source (Rohde & Schwarz SGS100A) with an internal IQ mixer to drive the readout resonators.We use a single-sideband mixing scheme to address multiple resonators coupled to the same feedline using a single microwave local oscillator (LO).The intermediate frequency (IF) signals are generated using two arbitrary waveform generator (AWG) channels (Keysight M3202) and are interfered with the LO using the IQ mixer to obtain an RF signal with multiple frequency components, each of which can address a different resonator.In order to avoid saturating the IQ mixer, we attenuate the output signal of the AWG channels using 10 dB attenuators.
At room temperature, the resonator response signal is amplified by a HEMT amplifier (MITEQ AMF-5D-00101200-23-10P) and sent through a band-pass filter (K&L 033F8) to filter out the TWPA pump tone.The signal is then downconverted using a heterodyne demodulation scheme [4] using an IQ mixer (Marki IQ-4509LXP) and the readout probe LO.The demodulated signals resulting from the interference between the readout output signal and the LO are then filtered using low-pass filters (Mini-circuits VLFX-300+), passing only the terms at different frequencies corresponding to the signal from each resonator.The two output signals are then digitized on two separate digitizer channels (Keysight M3102), and the signal for each frequency component is processed using a built-in field-programmable gate array (FPGA) to obtain the static in-phase, I, and quadrature, Q, signals.
We use the secondary LO output of the microwave source for the demodulation LO.In order to reduce the signal noise, we use two high-pass filters (Mini-Circuits VHF-6010+) in series.We use a combination of attenuators and an amplifier (Mini-Circuits ZVE-8G+) to ensure that the LO power input to the IQ mixer is in the 10-13 dBm range required for optimal mixer performance.
A 24-channel DC voltage source (QDevil QDAC-I) with in-line resistors generates the current required to apply a static flux to each qubit in our setup.In addition to this static biasing of the qubit, we can apply a baseband, fast-flux pulse to each qubit individually with an AWG channel (Keysight M3202), which allows us to tune the qubit frequency of the qubits on the timescale of nanoseconds.We attenuate the flux pulse produced by the AWG using a 10 dB attenuator in order to reduce the flux noise experienced by the qubit.
For individual X, Y control of each qubit, we use a microwave source (Rohde & Schwarz SGS100A) with an internal IQ mixer.We calibrate the I and Q quadrature offsets, the gain imbalance, and the skew of our IQ mixer by minimizing both the LO leakage and the signal level of the unwanted sideband.We use single-sideband mixing to generate our drive pulses and use an additional AWG channel to gate the RF output in order to mitigate the adverse effects of mixer LO leakage.In our setup, we use four different AWG channels to fully control each qubit: 1 channel for fast Z control, and 3 channels for X, Y control.
We combine the DC bias (used for tuning the static frequency of the qubit), baseband pulse (used for fast tuning of the qubit frequency), and RF pulse (used for charge driving the qubit) at RT.The DC and baseband flux signals are combined using a resistive combiner (Mini-circuits ZFRSC-42B-S+).In order to match the impedance of the resistive combiner to 50 Ω to prevent distortions in the baseband pulse, we use a homemade resistor box on the DC side, consisting of a 1 kΩ resistor in parallel with a 50 Ω resistor in series with a 10 µF capacitor connected to ground.The output flux signal is sent through a 300 MHz low-pass filter (Mini-circuit VLFX-300+) and is combined with the RF signal using a diplexer (QMC-CRYODPLX-0218).We use a 1 GHz low-pass filter (VLFG-1000+) with a flat response of approximately 40 dB of rejection throughout the qubit operation frequency range.Using this filter, we are able to filter out high-frequency sources of noise while maintaining the ability to apply high-fidelity single-qubit gates.
The advantage of this approach is that we only need a single microwave coaxial line for full control of each qubit.In addition, by combining DC and baseband flux signals at RT we are able to more efficiently characterize and compensate for the flux transients arising from signal combination with a fast oscilloscope or a digitizer.However, the impedance mismatch caused by the reflective low-pass filter at the mixing chamber gives rise to standing modes between the filter and the chip ground at the end of the flux line.As a result, the effective qubit-drive coupling depends on the drive frequency.We characterize this dependence by biasing the qubit at different frequencies and Rabi driving on resonance.
The active components used in our experimental setup are summarized in Table SI.

S2. SAMPLE
The 4 × 4 superconducting transmon array is fabricated using a 3D integration flip-chip process [5] (Fig. 1f in the main text).The qubits are located on a 4 × 4 mm qubit tier as shown in Fig. 1g, and the readout and control lines are located on a separate 5 × 5 mm interposer tier as shown in Fig. 1h.The two layers are separated using 3 µm superconducting indium bump bonds and silicon hard stops.
The interposer tier contains the flux bias lines and readout resonators that are respectively inductively and capacitively coupled to the qubits across the gap separating the two chips.The ground plane of the interposer layer is etched away in the regions corresponding to the location of each qubit, with an additional 10 µm gap to avoid unwanted capacitance between the qubits and the interposer ground plane.
The qubit tier comprises a 2D lattice of capacitively-coupled single-ended transmon qubits [6] arranged in 4×4 array.The capacitive coupling between each transmon is facilitated using an intermediary piece of superconductor, resulting in a nearest-neighbor qubit-qubit exchange interaction with average strength J/2π = (5.89± 0.4) MHz, measured at qubit frequencies of 4.5 GHz.The ground plane of the qubit tier is etched away in the regions corresponding to the location of the signal line on the interposer tier, with an additional gap equal to the width of each line.
Multiplexed, dispersive qubit-state readout is performed through individual capacitively coupled coplanar resonators.The sample contains two readout feedlines, each containing a single-port Purcell filter coupled to eight readout resonators.We use λ/4 readout resonators for all qubits along the edge of the lattice to minimize their footprint on the sample.The central qubits are read out with a λ/2 resonator.This is done to minimize the parasitic couplings; the central qubits' resonator crossing is located in the middle of the λ/2 resonator at its voltage node.The Purcell filter is designed to have a large bandwidth of 0.5 GHz at a resonance frequency of 6.3 GHz, such that the readout resonators can be distributed over a ∼ 400 MHz band while limiting Purcell decay to a rate ≲ 1/(300 µs).The resonators have a measured average linewidth of κ/2π = (1.17± 0.3)MHz, with an average dispersive shift of χ/2π = (0.87 ± 0.11)MHz characterized with each qubit biased at 4.5 GHz.
The 3D integration process allows us to route each flux bias line directly above the corresponding superconducting quantum interference device (SQUID).This positioning allows us to attain an average mutual inductance between the qubit SQUID loops and their respective flux bias lines of 1.15 pH while reducing the SQUID area to 4 × 4 µm compared to planar transmon arrays [7,8].Decreasing the area of the SQUID reduces the susceptibility of the qubit to flux noise, the offset due to uncontrollable magnetic fields in the environment, and crosstalk.We measure more than an order of magnitude reduction in the crosstalk levels of the flip-chip sample compared to a planar 3 × 3 array of transmons [7,8].
The sample is fabricated on a silicon substrate by dry etching an MBE-grown, 250-nm-thick aluminum film in an optical lithography process, forming all larger circuit elements such as the qubit capacitor pads, resonators, and the signal lines for qubit readout and control.The qubit SQUID loops are fabricated with an electron beam lithography process and a double-angle shadow evaporation technique [9] to form the Josephson junctions.We use a wire width of 5 µm for the SQUIDs to minimize flux noise from local magnetic spin defects on the SQUID surface and interfaces [10].
Detailed sample parameters for individually isolated qubits are summarized in Table S2.We show the maximum transmon transition frequencies ω max q at the upper flux insensitive point, the qubit anharmonicities U (measured at ω max q ), the readout resonator frequencies ω res , the probabilities F ij of measuring the qubit in state i after preparing it in state j, the readout assignment fidelity (F gg + F ee )/2, the measured T 1 s and T * 2 s at qubit frequencies around the bias point used in the experiment.We anticipate that the reported readout assignment fidelities are limited by the state preparation error due to the initial thermal population of the qubits.The relatively low F gg is a result of the high thermal population caused by noticeable heating at the MXC stage when DC biasing the qubits.This heating is caused by sending a DC current through the microwave cable between the still and MXC stages, dissipating approximately ∼ 500 µW of power at the MXC when all the qubits are biased at 4.5 GHz.In addition, we report the individual and simultaneous single-qubit randomized benchmarking (RB) fidelity for the gates used for tomographic readout.The microwave charge pulses for the gates are applied via the flux lines.Due to the significant RF crosstalk between the drives targeting QB5 and QB10, we do not calibrate and apply simultaneous gates for tomographic readout on these qubits.

S3. CONTROL AND READOUT CALIBRATION A. Flux crosstalk calibration
Flux-tunable transmon qubits comprise two Josephson junctions in a SQUID loop in parallel with a shunting capacitor.We can control the effective Josephson energy E J of the SQUID, and therefore the transition frequency between the ground and the first-excited state of the transmon, by threading magnetic flux through the SQUID.The frequency of the transmon in response to an applied external flux Φ ext is approximately given by [6]: where E C is the transmon charging energy and )/h is the maximum qubit frequency.The asymmetry parameter d of the SQUID junctions is given by d = |(E J,2 − E J,1 )/(E J,2 + E J,1 )|, where E J,1 and E J,2 are the Josephson energies of the two SQUID junctions.The above equation is an approximation because it holds only in the limit at the upper flux-insensitive point, the qubit anharmonicities U (measured at ω max q ), the readout resonator frequencies ωres, the probabilities Fij of measuring the qubit in state i after preparing it in state j, the readout assignment fidelity (Fgg + Fee)/2, the measured T1s, T * 2 s.The reported readout assignment fidelities are limited by the state preparation error due to the initial thermal population of the qubits.In addition, we report the individual and simultaneous single-qubit randomized benchmarking (RB) fidelities.Due to the significant RF crosstalk between the drive targeting QB5 and QB10, we do not apply single-qubit gates to these qubits in our experiments.
We thread flux through the SQUID by applying a current through a local flux line that terminates near the SQUID.We use a voltage source applied over a series resistor of resistance 1 kΩ at room temperature to generate a stiff current source for each flux line.There is a linear relation between the voltage applied to a flux line and the magnetic flux experienced by the SQUID: where V Φ0 is the voltage required to tune the qubit by one magnetic flux quantum Φ 0 , and Φ offset is a flux offset due to non-controllable sources of static magnetic field.
For an array of transmons, the applied current in a flux line targeting one SQUID may thread unwanted magnetic flux through other SQUIDs.We can model this flux crosstalk as a linear process, where the fluxes ⃗ Φ ext experienced by the SQUIDs are related to the voltages ⃗ V applied to the flux lines: where V Φ0 is a diagonal matrix with V Φ0 i,i corresponding to the V Φ0 of qubit i. S is the flux crosstalk matrix, with S i,j = ∂V i /∂V j representing the voltage response of qubit i to a voltage signal applied to qubit j.In this representation, the diagonal elements of S are 1.By characterizing S, we can compensate for the crosstalk and set the qubit frequencies precisely.Response of Qubit i 100 0.6 -0.1 -0.4 0.9 1.0 -0.7 -0.8 0.2 0.0 -0.8 -0.7 -0.4 -0.5 -0.6 -0.7 1.6 100 0.2 -0.2 0.8 1.2 -0.8 -0.9 0.4 0.1 -0.9 -0.9 -0.3 -0.4 -0.7 -0.9 -0.5 -0.1 100 1.6 -0.9 -0.7 1.1 0.8 -0.8 -0.9 0.2 0.4 -0.7 -0.6 -0.5 -0.2 -1.5 -0.1 0.9 100 -0.8 -0.8 0. The first step of calibrating flux crosstalk is characterizing the transmon spectrum of each qubit in the array.We do this by sweeping the voltage applied to the flux line and measuring the qubit frequency.By fitting this data with Eq. 1 and Eq. 2, we extract the transmon spectrum parameters f max , E C , and d, and the voltage to flux conversion parameters V /V Φ0 and Φ offset .
Next, we calibrate flux crosstalk with a learning-based protocol described in [11].We start with an initial guess crosstalk matrix, which can either be an estimate from a previous measurement or the identity matrix.
We select a quasi-random vector of target frequencies ⃗ f , subject to a few constraints.First, we require that the frequency for each qubit falls in the range of 100 MHz below its sweet spot to 1 GHz below its sweet spot.Second, to reduce frequency shifts due to resonant exchange interaction between qubits, we require that the frequency detuning between neighboring qubits is at least 200 MHz and that the frequency detuning between any two qubits in the array is at least 50 MHz.We use our initial S and equations Eq. 1 and Eq. 3 to solve for the voltages ⃗ V needed to target ⃗ f .We apply ⃗ V , then measure the qubit frequencies via spectroscopy.Even though we detune our qubits, they still experience frequency shifts.We use previously characterized qubit-qubit couplings to calculate the uncoupled qubit frequencies.Finally, we convert the uncoupled frequencies to fluxes experienced.
By repeating this process M times, we obtain a size-M training set of applied voltages and measured fluxes: Φ meas,i } i=1:M .We learn our crosstalk matrix using the crosstalk relation (Eq.3), by leveraging the fact that the estimated flux (V Φ0 ) −1 S ⃗ V i + ⃗ Φ offset ought to equal the measured flux ⃗ Φ meas,i .We learn S row-by-row.We learn the elements of the k th row of S by minimizing the mean-squared-error cost function in a gradient descent optimizer (we use the L-BFGS optimizer in PyTorch).The optimization will converge as the estimated fluxes approach the measured fluxes.Once we have repeated this minimization for all rows, we have our optimized S. In Fig. S2a, we report the DC flux crosstalk matrix for the 16-qubit device.This same calibration procedure is utilized to calibrate crosstalk for fast flux pulse control, as well.In order to learn the fast flux crosstalk for each qubit, we tune that qubit to a target frequency using a 100 ns fast flux pulse and measure the qubit frequency via spectroscopy.Then, using the same pulse amplitude for the target qubit, we apply flux pulses with random amplitudes through the other flux lines and measure the change in the target qubit's frequency.We use the training set consisting of the voltages applied and the changes in the target qubit frequency for learning the fast flux crosstalk matrix.In Fig. S2b, we report the fast flux crosstalk matrix for the 16-qubit device.

B. Time-domain flux pulse shaping and transient calibration
Baseband-frequency flux pulses experience distortions as they travel to the sample.The distortion caused by roomtemperature components can be characterized and compensated for using a fast oscilloscope with a sampling rate better than 500 MHz.The distortion introduced by the components in the cryostat is generally temperature-dependent and needs to be characterized while the fridge is operating.By using the qubit as a sensor, we can characterize the distortion of a square flux pulse (Z pulse) using Ramsey-type measurements [12].The step response for the signal generated by the AWG, s AWG (t), the signal reaching the qubit, s qubit (t), and the time-dependent distortion, n(t), are related by To characterize n(t), the target qubit is biased to its sweet spot using static flux and excited with a π 2 -pulse around the x-axis followed by a Z pulse for time t.The dynamic frequency change of the qubit as a response to the Z pulse can be extracted by measuring both ⟨X⟩ and ⟨Y ⟩, denoting the expectation values of the qubit state projected along the x-axis and y-axis of the Bloch sphere.The phase accumulated by the qubit during the Z pulse can be described by ϕ(t) = t 0 ∆(t ′ )dt ′ , where ∆(t) = ∆ 0 + δ(t) is the frequency detuning, consisting of the target detuning ∆ 0 and the added frequency distortion δ(t).The ⟨X⟩ and ⟨Y ⟩ measurements are used to obtain cos ϕ(t) and sin ϕ(t), respectively, which we use for constructing the phasor e iϕ(t) = cos ϕ(t) + i sin ϕ(t).We can then extract the frequency distortion as δ(t) = d dt arg(e iϕ(t) ).By using the qubit spectrum, we map δ(t) to the distortion in the signal reaching the qubit n(t).
In order to obtain an analytical expression, we fit the extracted n(t) with a sum of typically three or more exponential damping terms with time constants τ i and settling amplitudes A i : In our experiments, we use the analytical expression of n(t) to pre-distort the output signal of the AWG in order to compensate for the system transients.
One method of analytically pre-distorting the signal is using Fourier transformations.Using the signal step-response h(t) = Θ(t)(1 + n(t)), we can obtain the impulse-response of our system: We can solve for the pre-distorted signal needed to make the qubit feel a constant flux using the convolution relation s qubit (t) = s AWG (t) * I(t), therefore, where We can calculate an expression for F{I(t)}: In practice, in order to accurately pre-distort the pulse shape using Fourier transformations we need to zero-pad the waveform such that the total length of the waveform is at least 5 × max i (τ i ).As a result, when correcting for long-timescale transients with τ ∝ O(100 µs) this procedure becomes inefficient.
Instead, we use an infinite-impulse-response (IIR) filter for pulse pre-distortion [12] with the difference equation: The timing for the Z pulse applied through each line is calibrated by sweeping the delay time of a Z pulse by tuning the qubit to the sweet spot while applying an XY pulse through the common resonator feedline at that frequency.(b) The timing for the XY pulse applied through each local line is similarly calibrated by applying a Z pulse with a constant delay characterized in (a) and sweeping the delay time of the XY pulse.
where y[n] is the output signal, x[n] is the input signal, N is the feed-forward filter order, and M is the feed-back filter order.A first-order IIR filter (N = M = 1) can be used to correct the exponential transients.We find the IIR coefficients a 0 , a 1 , b 0 , and b 1 for a single exponential pole with amplitude A and time-constant τ : where α = e −1/fsτ depends on the AWG sampling rate f s .We use the SciPy library in Python to apply the IIR filter with the coefficient calculated above to the waveforms in our pulse generation software.For systems with more than one exponential transient term, we concatenate multiple IIR filters, each corresponding to a term.

C. Pulse alignment
Differences in the pulse path lengths and the delays in control electronics cause discrepancies in the time that flux (Z) and charge (XY ) pulses reach the target qubit.To accurately execute the desired pulse sequence, we characterize the relative delay of each signal line and adjust pulse timings in software to compensate.
We first characterize the delay of a Z pulse applied to each qubit with respect to a common XY drive through the resonator feedline.An XY pulse through the resonator feedline reaches each qubit roughly simultaneously; it, therefore, serves as a common reference for pulse delay timing.To calibrate the Z pulse delay of a particular qubit, we first isolate the qubit from the rest of the lattice and away from the flux sweet spot.We then apply a short 13 ns Z pulse with the amplitude required to tune the qubit to the sweet spot, and an XY pulse of the same length at the sweet-spot frequency.The amplitude of the XY pulse is less than that required for a π-pulse.We measure the qubit resonator response as a function of the Z pulse delay (Fig. S3a); a peak in the measured response indicates the delay at which the qubit is maximally excited, which occurs when the two pulses are aligned.(When the two pulses are misaligned, part, or all, of the XY drive is off-resonant with the qubit frequency, resulting in a lower qubit excited state population.)By fitting the data in Fig. S3a to a Gaussian curve, we find the optimal Z delay of 17 ns with respect to the XY pulse applied through the common line.Next, we calibrate the delay for XY pulses sent to each qubit through its respective local line, using the calibrated Z pulse delays for reference.Similar to the procedure described above, we tune the qubit to the sweet spot using a Z-pulse and apply an XY pulse, now through the local line, at the qubit sweet spot frequency.Now with the Z pulse timing fixed at the calibrated value, we sweep the delay of the local XY pulse and extract the optimal local XY pulse delay with respect to the common XY pulse timing, as shown in Fig. S3b.

D. Readout optimization
We use frequency multiplexed heterodyne mixing to read out multiple resonators coupled to the same feedline with a common local oscillator (LO) frequency [4].To maximize the readout fidelity, we maximize the separation in in-phase/quadrature (I/Q) space between the readout signal conditioned on the qubit ground and excited state by optimizing the frequency and amplitude of each sideband tone.
In order to distinguish between the single-shot clusters resulting from different qubit states we first take single-shot measurements with the qubit prepared in the |0⟩ state and the |1⟩ state by applying a π−pulse to the qubit.Using this data, we train a support vector machine (SVM) [13] to find the hyperplane that classifies each single-shot point as 0 or 1.The hyperplane parameters are calculated in a manner that maximizes the probability of correctly classifying each data point.We quantify the readout quality using the readout fidelity: where f 00 = P (0| |0⟩) is the probability of correctly classifying the qubit |0⟩ state as 0, and f 11 = P (1| |1⟩) is the probability of classifying the |1⟩ state as 1 after readout.Fig. S4a illustrates an example set of single-shot data acquired by preparing the qubit in the |0⟩ state (blue) and in the |1⟩ state (red) and the corresponding line used to discriminate between the two states.For this dataset, the readout fidelity is F RO = 98.0%, with f 00 = 99.6% and f 11 = 96.5%.The discrepancy between f 11 and f 00 is a result of qubit relaxation during the readout process.
In addition to using F RO to quantify the readout quality, we also consider cluster separation.We use a Gaussian Mixture Model (GMM) probability distribution to group the single-shot data into two clusters, independent of the intended state.We fit the mean (µ) and the covariance matrix Σ for each Gaussian distribution representing a cluster using the expectation-maximization (EM) algorithm [13].A multivariate Gaussian distribution in the I/Q plane with mean µ and covariance matrix Σ has density: In our model, we use the same covariance matrix Σ for both clusters.We use the spatial overlap between the Gaussian distributions to quantify the cluster separation by defining the separation fidelity F sep as In Fig. S4b, we group the readout single-shot readout into two clusters.The ellipse around each cluster encloses the area corresponding to two standard deviations away from the center of the cluster.Using Eq. 11 we obtain a separation fidelity of F sep = 99.9% between the two readout clusters.
To optimize the readout in our calibration procedure, we vary the sideband modulation pulse frequency and amplitude used for multiplexed readout.Both the readout fidelity F RO and the separation fidelity F sep can be used as the optimization metric.We measure the readout optimization landscape by sweeping the readout frequency offset and readout amplitude for a given qubit and calculating 1 − F RO (Fig. S5a) and 1 − F sep (Fig. S5b).The optimization landscape using both metrics is convex, which allows for the use of a conventional optimizer, such as Nelder-Mead, to find the optimal readout sideband frequency and amplitude efficiently.The landscape of the cluster separation fidelity is more sensitive to changes in the readout parameters compared to the readout fidelity, which is a result of the limit imposed on the readout fidelity caused by qubit relaxation during readout.Therefore, separation fidelity is a better metric to use for the optimal readout sideband frequency and amplitude while keeping the readout time fixed.It is important to limit the maximum readout amplitude to ensure that the readout remains in the dispersive regime In the dispersive readout scheme, with the readout resonator frequency above the transmon frequency, we find the magnitude of the optimal readout frequency offset to be roughly equal to the measured qubit-resonator dispersive shift χ as predicted for χ ≲ κ/2 [4].
Using these optimal parameters, we vary the readout integration time to maximize the readout fidelity (Fig. S5c).Increasing the integration time results in an increase in the single-shot cluster separation, however, the readout fidelity suffers more from qubit relaxation.Therefore, the optimal integration time is selected in a way to balance these two effects by maximizing F RO .

E. State preparation and measurement (SPAM) error mitigation
In order to mitigate SPAM errors from our measured observable, we use a stochastic β-matrix (often also referred to as a T matrix in the literature) for each qubit, where where f 00 and f 11 are the probabilities of correctly preparing and measuring the qubit in the |0⟩ and |1⟩ states respectively.There are two main sources of our experimental error captured by this matrix: incoherent thermal population of the qubit prior to the experiment sequence, and readout errors.
In the absence of readout crosstalk, the β-matrix for an n-qubit system is the tensor product of the β-matrix for each individual qubit: β = β 1 ⊗ ... ⊗ β n .By inverting this matrix, and applying it to the measured bitstring probabilities we can correct for SPAM errors in the measured observable [14].

F. Idling frequency layout
During the emulation of the Bose-Hubbard Hamiltonian, all qubits are brought on resonance to allow particle exchange.During state preparation and readout steps of the experiment sequence, the qubits are detuned from one another to turn off particle exchange.The capability to apply simultaneous high-fidelity single-qubit gates to qubits at their respective idling frequencies is necessary for performing tomographic readout, which is required to measure correlations and state purities.The idling frequencies of the qubits during the state preparation and readout steps of the experimental sequence satisfy the following criteria: • Each qubit is detuned from its nearest neighbors by at least 300 MHz to prevent particle exchange during readout.
• Each qubit is detuned from its next-nearest neighbors by at least 40 MHz.
• Qubit frequencies are selected to minimize charge-driving crosstalk effects to enable high-fidelity simultaneous single-qubit gates.
• Qubit 12 is biased at 4.5 GHz (qubit 12 suffered from elevated pulse distortion due to a partial break in its flux line).
• Each qubit frequency is detuned from any two-level system (TLS) defects.
The idling qubit frequency layout used in our experiment is shown in Fig. S6.At these frequencies, we apply singlequbit gates on all qubits, excluding qubits 5 and 10, for which the qubit was too weakly coupled to the line used for charge driving.The randomized benchmarking (RB) fidelities for the calibrated gates are reported in Tab.S2, with an average individual RB fidelity of 99.8% and an average simultaneous RB fidelity of 99.6%.

S4. HAMILTONIAN CHARACTERIZATION
The Hamiltonian of the driven system is described by as shown in Eq. 5 of the main text.The particle exchange strengths J ij and the amplitudes and phases of the relative drive couplings α i are determined by the circuit layout, and remain constant in our experiments.To accurately simulate the behavior of our lattice, we carefully characterize each of these parameters.

A. Particle exchange strengths
We measure the particle exchange strengths J ij for nearest and next-nearest neighboring qubits in the lattice (that is, for qubits i and j of relative Manhattan distance 1 and 2).
Each value J ij is measured via the resonant iSWAP rate between qubits i and j at ω com /2π = 4.5 GHz.With all other qubits detuned, we first add an excitation to qubit i by applying a π-pulse.Qubits i and j are then brought on resonance at frequency ω com .After allowing the qubits to interact for time t, the qubits are once again detuned and measured in the z basis.The coupling strength J ij is determined by the time T required for a full excitation swap, J = π/2T [15].
A histogram of the measured values J ij for (next-)nearest neighbors are shown in Fig. S7(a(b)).The average nearest coupling in our device is J/2π = 5.9(4) MHz at ω com /2π = 4.5 GHz.The next-nearest neighbor couplings are typically in the range of J/10 to J/40.The qubits at the center of the lattice exhibit a larger coupling to their next-nearest neighbor, which can be attributed to larger stray capacitances.

B. Drive coupling amplitude
To determine the coupling amplitude g i = Ω × |α i | between qubit i and the common drive, we use a Rabi measurement.We drive each qubit on resonance at ω com /2π = 4.5 GHz, with all other qubits far-detuned to avoid particle exchange.We sweep the drive amplitude A and the drive pulse duration, and measure the Rabi rate at each drive amplitude (Fig. S8a).Next, we find the Rabi rate for different drive amplitudes and use a linear fit to extract g i /A (Fig. S8b).We repeat this procedure for all qubits in our lattice to find the Rabi rate dependence on the common drive amplitude.By asserting that |α| = 1 for the qubit with the strongest coupling to the common drive (qubit 11), we report the relative drive coupling magnitudes for all the lattice qubits in Fig. S8c (see Tab. S4 C for the values).The common drive coupling strength Ω is therefore related to the drive amplitude A through the linear relationship We notice that when driving all 16 qubits simultaneously while on resonance, we need to apply a +7.5% correction to the relation above in order for all of our numerical simulations to match the experimental data.

C. Drive coupling phase
We use a two-step approach to characterize the coupling phase ϕ i = arg α i between the common drive and each qubit.We need to calculate 15 coupling phases ϕ i relative to ϕ 1 (the global phase is irrelevant).Single-qubit experiments cannot provide information on coupling phases, so to determine them, we need to conduct drive coupling characterization experiments on multi-qubit subsystems of the entire lattice.
We initially characterize the relative drive coupling phase between different pairs of neighboring qubits using an interferometric approach.For such experiments, we isolate a qubit pair by far-detuning the rest of the lattice.We begin by biasing qubit "A" at ω com /2π = 4.5 GHz and applying a π 2 -pulse through the common line while detuning qubit "B" to prevent them particle exchange.After this step, the two-qubit system is in the state where ϕ A is the phase of the common drive coupled to qubit "A".Afterward, we bring qubit "B" on resonance with qubit "A", allowing the two qubits to undergo an iSWAP interaction for time t.After a full iSWAP period, the state of the system is Next, we detune qubit "A", while keeping qubit "B" at ω com and apply a π 2 -pulse through the common line to qubit "B" resulting in the state where ϕ B is the phase of the common drive coupled to qubit "B".By sweeping the phase ∆ϕ of the drive and measuring the population on qubit "B" we can find ϕ A − ϕ B .The excited state population of the qubit will depend on ∆ϕ Hence, when ∆ϕ = (ϕ A − ϕ B ) the measured population is p(|1⟩) = 1/2 with a negative slope.We show the pulse sequence for this measurement in Fig. S9a, and representative measurement results using qubits 9 and 10 in Fig. S9b.The gray dashed line indicates ∆ϕ = (ϕ 9 − ϕ 10 ).Using these relative phases measured between the different qubit pairs, we calculate an absolute phase for each qubit with reference to the drive phase on qubit 1, with ϕ 1 = 0.
We then take a second step to improve the accuracy of characterizing the drive coupling phases.We isolate different 2 × 2 plaquettes within our lattice and set the corresponding qubit frequencies to ω com .We apply a weak common drive with strength Ω = J/2 and measure the population on each qubit within the plaquettes as a function of time.The phase difference between the adjacent qubits extracted from the respective 2 × 2 plaquette fits.The two numbers in grey were unused because those bonds were captured through a higher-quality fit in the neighboring plaquette.The two red-shaded plaquettes were unused because of poor fit quality.
By comparing the measured qubit populations to numerical simulations, we fine-tune the drive coupling phases.To accomplish this, we define the cost function where P data i (t) is the measured population of qubit i after time t and P sim i (t, ⃗ ϕ) is the corresponding simulated value, as parameterized by the three relative drive coupling phases ⃗ ϕ in the plaquette.Starting with the characterized drive coupling phases using the pairwise interferometry experiments, we use a gradient-descent-based optimizer to find phase values that minimize the cost function.
We utilize the JAX [16] and OPTAX [17] packages in Python to perform simulations and optimization.We leverage automatic differentiation to compute gradients and determine parameter updates with an Adam optimizer.The initial values provided by the first drive coupling phase calibration step are crucial for gradient-descent-based minimization of the cost function since C is non-convex in ⃗ ϕ.We report the experimental data from the different plaquettes and the numerical simulations with the optimized drive coupling phases in Fig. S9c.Generally, we find an excellent agreement between the data and numerical simulation.However, for two 2 × 2 plaquettes, formed by qubits (5, 6, 9, 10) and qubits (7,8,11,12), we observed poor convergence of the optimizer, hence, we exclude the results from those two plaquettes from estimating the phases.
The second step of characterization overdetermines the relative phases: we performed the second step for seven 2×2 plaquettes.Each such measurement yields three independent phase parameters, producing a total of 21 extracted relative phases.Yet there are only 15 physically independent phase values.In order to estimate the most likely drive coupling phases, we take the quality of each fit for each plaquette into consideration (see Fig. S9d).In Tab.S4 C we report the final characterized values for ϕ i = arg(α i ) that we use for the numerical simulations of our 4 × 4 lattice.
To emphasize the importance of accurately characterizing the drive coupling phase to numerically simulate the behavior of our driven lattice, we consider the time dynamics of the average number of excitations ⟨n⟩ in the lattice under a resonant drive.In Fig. S10 we compare experimental data for drive strength Ω = J/2 with numerical simulations with (i) all phases set to zero, (ii) using phases characterized with the pairwise interferometry experiments and (iii) with fine-tuned phase values.We observe that with the fine-tuned phases the numerical simulations of our 4 × 4 lattice are in excellent agreement with our experimental results.

S5. HARD-CORE BOSE-HUBBARD MODEL ENERGY SPECTRUM
We study the energy spectrum of the hard-core Bose-Hubbard Hamiltonian ĤHCBH with uniform site energies (ϵ i = 0) described in Eq. 4 of the main text.We use the measured qubit-qubit coupling strengths, including nearest neighbor and next-nearest neighbor exchange interactions, characterized at 4.5 GHz.In order to efficiently diagonalize the Hamiltonian of our 16-qubit lattice, we first project the 2 16 × 2 16 matrix into a subspace spanned by a fixed particle number n.This projection can be accomplished using a 2 16 × 16 n matrix U , where the columns of U are all the permutations of the vectors corresponding to the product states with exactly n particles.The Hamiltonian corresponding to the n particle subspace is therefore calculated through where Ĥ′ n is an 16 n × 16 n matrix.By diagonalizing Ĥ′ n we obtain the eigenenergy distribution for different particle numbers n as shown in Fig. S11a.The largest reduced Hamiltonian is Ĥ′ n=8 (describing the subspace of n = 8) and its width is 2 16 / 16  8 ≈ 5 times less than the full Hamiltonian.Therefore, we expect over 100× speed-up when diagonalizing the Hamiltonian of each particle number subspace independently compared to diagonalizing the full system Hamiltonian.We notice a mild skew towards the higher positive energy eigenenergies in the energy spectrum of our lattice obtained numerically (see Fig. S11a).In contrast, the eigenenergies of the system Hamiltonian excluding the NNN coupling terms are symmetric around zero energy (Fig. S11b).In Fig. S11c we show the difference between the magnitude of the highest eigenenergy E max and the lowest eigenenergy E min for the different particle-number subspaces.We notice that for all particle-number subspaces, the eigenenergy distribution is skewed in the positive direction for our system Hamiltonian, which includes NNN couplings, whereas in the absence of NNN coupling, the eigenenergy distribution is symmetric around E = 0. Therefore, we conclude that the NNN coupling present in our lattice causes a skew in the eigenenergies of the system.For each constant-particle-number subspace of the HCBH Hamiltonian, we observe a variation in the geometric entanglement from the edge to the center of the spectrum [18].To illustrate this variation, we report the average subsystem entropy as a function of volume for states at the edge and the center of the energy band of subspaces with n = 5, 6, 7, 8 particles in Fig. S12a.The states at the center of the energy band exhibit a distinct Page curve, while the entropy of the states at the edge of the energy band shows a weak dependence on volume.Furthermore, in Fig. S12b we show the geometric entanglement ratio s V /s A and notice the same trend between the states at the center and edge of the energy band for the subspaces designated by the different number of particles.The geometric entanglement behavior is consistent across different particle-number subspaces, allowing us to probe the entanglement scaling across the many-body spectrum using a superposition of different eigenstates.
Instead of preparing specific eigenstates, we can use a superposition of eigenstates in different regions of the energy spectrum to probe the entanglement properties in a many-body system.We accomplish this by applying a weakly driving a uniform lattice with some detuning δ from the frequency of the qubits.
The hard-core Bose-Hubbard Hamiltonian for a lattice with N sites can be represented in the eigenstate basis as where |n, ϵ⟩ is the eigenstate and ρ(n, ϵ) is the density of states with n particles and energy ϵ.In this basis, the driving operator Σ can be represented as For a weak drive, the driving operator will couple eigenstates that are separated in energy by detuning ϵ ′ − ϵ = δ.Therefore, we can approximate the operators as a combination of energy-raising operators with where ϵ n = ϵ × δ.Each Âϵ couples states on the Hamiltonian spectrum shown in Fig. S11a that are on a line with slope δ and intersection n = 0 at energy ϵ.For a lattice initialized with no particles, the only relevant operators are Â0 and Â † 0 .Therefore, for times t much shorter than 1/Ω, the state of the system is approximately which is similar to a multi-mode coherent state.Due to exchange and on-site interactions, the analogy to coherent states becomes weaker at longer times as the number of photons becomes larger.In this section, we discuss the impact of the characterized drive coupling phases on our experiments using numerical simulations.First, we consider the impact of the drive coupling phase on the formation of entanglement entropy within a driven lattice.In Fig. S13a we show the simulated time evolution of the average entanglement entropy of the 8-qubit subsystems within the 4 × 4 lattice when driven on-resonance with drive strength Ω = J/2.We see that the entropy within the lattice driven with the characterized phases reaches equilibrium much faster than a driven lattice with uniform phases.
Next, we study the correlation lengths within the states prepared with and without the characterized drive phases (Fig. S13b).To ensure that the prepared states are in equilibrium, we simulate evolution of the lattice for t = 10/J under the drive with the characterized coupling phase values, and t = 50/J for the drive with uniform phases.We observe that coherence lengths follow the same trend in both preparation scenarios, with some minor deviations.Finally, we simulate the effect of the drive phase on the correlation matrix in Fig. S13c.The correlation matrices exhibit the same general pattern, however, the correlations between neighboring sites are stronger for states that are prepared using the drive with the characterized coupling phases.In order to quantify the total entanglement formed, we use the normalized average purity of all the qubits in the lattice E gl = 2 − 2 N i Tr(ρ 2 i ), which serves as a global entanglement metric [19].This quantity can be measured with low experimental overhead.We measure E gl by reconstructing the density matrix for each qubit in the lattice via single-qubit tomography.In Fig. S14a we report the measured E gl after driving the lattice for time t with a drive detuning δ.In Fig. S14b, the close agreement between experimental data and numerical simulations, which do not take into account decoherence, suggests that the entanglement formed is between different sites of the lattice, rather than with the uncontrolled environmental degrees of freedom related to decoherence.
While we use E gl to probe the formation of entanglement within our lattice, we cannot distinguish the nature of the entanglement (i.e.whether it is short-ranged or long-ranged).The dynamics of the entanglement formation in the driven lattice reach a steady state after t ≈ 10/J.We note that as the drive detuning δ gets larger, the steady state value for E gl decreases.This decrease can be attributed to a reduction in the average number of excitations in the lattice with an increase in δ (see Fig. S28).

S9. TOMOGRAPHY SUBSYSTEMS FOR STUDYING THE ENTANGLEMENT BEHAVIOR
In order to tomographically reconstruct the density matrix for a 6-qubit subsystem, we need to measure 3 6 = 729 Pauli-strings (e.g.ZXXY ZZ).We use a maximum-likelihood estimation (MLE) routine implemented in Qiskit [20] to reconstruct the density matrix using the measured Pauli-strings.With the aid of our high-fidelity simultaneous single-qubit gates and readout, we can efficiently measure the Pauli-strings needed to reconstruct the density matrix for multiple 6-qubit subsystems using a single round of measurements.Our approach involves initially listing the complete collection of Pauli-strings needed to reconstruct the 6-qubit subsystem (3,4,7,8,11,12).We then pair up each of the qubits in the rest of the lattice with a qubit in the subsystem so that their measurement basis is consistent for every measurement.To illustrate how we match the measurement basis of our qubits, we utilize Fig. S15a, in which we use matching colors to indicate the qubits that have identical Pauli operators in each Pauli-string.Qubits 5 and 10 are excluded from tomographic measurements as applying single-qubit gates to these two qubits causes substantial drive crosstalk affecting the other qubits in the system.
Next, we identify all the subsystems that can be reconstructed using the measured Pauli-string.We include all subsystems of size one through six that are nearest-neighbor connected and in which each qubit possesses a distinct color (according to Fig. S15a), yielding 163 unique subsystems.We report a list of all these subsystems, sorted by size V in Table S9.In Fig. S15b we report the distribution of the area (A X ) and the volume (V X ) of the subsystems used for studying the entanglement scaling behavior in the main text.In Fig. S16, we provide a visualization of the entropy extracted for all subsystems as a function of both area and volume.
In our experiments, we measure the expected value of each Pauli string with 2000 shots.We note that the number of samples we can extract for each Pauli string depends on the subsystem size.To see this, notice that measurements of the 6-qubit Pauli strings XY ZXY X, XY ZXY Y , and XY ZXY Z all include measurements of the 5-qubit Pauli string XY ZXY ; therefore, while measuring the former we get three times more measurements of the latter.In general, our measurements yield 2000 × 3 6−V samples of each volume V subsystem included in the coloring shown in Fig. S15a.

S11. MEASUREMENT SAMPLING STATISTICS
In the main text, we demonstrated excellent agreement in most cases between the entropy extracted from measurements and from simulations, yet we noted a slight deficit of the measured entropy for large subsystems in the volume-law regime.We attribute this deficit to be a consequence of sampling statistics during a finite number of measurements.In this section, we explain this effect, reproduce the effect in simulations that include measurement sampling statistics, and discuss mitigation strategies.Full-state tomography of a subsystem X containing V X sites involves measurement of Pauli strings i∈X σ αi i for all combinations of Pauli operators α i ∈ {x, y, z}.For each Pauli string, we aim to accurately determine the distribution of measurement outcomes, of which there are 2 V X .For larger subsystems, the number of possible measurement outcomes is large; and as the state being measured approaches infinite temperature (an ideal volume law), the distribution of measurement outcomes approaches a uniform distribution (since at infinite temperature all states are equally likely).In these limits, the number of measurements required to accurately sample the outcome distribution becomes large.
The area law states generated when |δ|/J is larger have lower absolute effective temperatures meaning they feature far-from-uniform distributions of measurement outcomes.Reconstruction of these states is therefore less sensitive to finite sampling statistics.This observation is commensurate with the results of Ref. [21], where only 3 × 10 3 samples per Paul string were sufficient to accurately reconstruct 10 qubit GHZ states (which have area-law entanglement scaling).
To quantify the impact of the number of samples n s on the extracted entropy, we take a Monte Carlo approach.Here, we consider the coherent-like state prepared at δ = 0 and Ω = J/2 (a volume-law state), and begin by obtaining the final state via a decoherence-free simulation on a classical computer.For each subsystem and each Pauli string, we then sample from the distribution of bitstring measurement outcomes n s times.We reconstruct the subsystem density matrices from these samples and compute their entropy S 2 .Density matrix reconstruction used the same maximum-likelihood estimation routine as was used to reconstruct density matrices from experimental data.
Results are shown in Fig. S18 for n s ranging from 50 up to 2 × 10 4 .Density matrix reconstruction without maximum-likelihood estimation is shown for comparison.For low n s , sampling bias causes a biased reconstruction of the distribution of measurement outcomes, resulting in a deficit of the extracted entropy.The extracted entropy increases and eventually saturates at the correct values as n s increases.The value of n s needed to accurately extract the subsystem entropy grows exponentially in subsystem volume.While n s = 2 × 10 3 was used for volume 6 subsystems in the present experiment, these simulations show that n s ≳ 10 4 is needed to accurately extract the entropy of volume-law states for subsystems of volume 6.
Results from Monte Carlo simulation of measurement sampling effects are compared to experimental data in     qubits.We observe that as the system size becomes larger, the contrast in s V /s A values extracted using different maximum subsystem sizes remains consistent.It is worth noting that as V max becomes larger the s V /s A (r=0.1) s V /s A (r=0.8) ratio decreases, however, for a given V max the value is consistent for various lattice sizes.This observation suggests that our technique involving measuring the entanglement entropy of subsystems of finite size in a lattice can faithfully distinguish between area-and volume-law states regardless of the overall size of the system.
The numerical simulations discussed in this section confirm that we can study the scaling of entanglement within larger lattices by performing tomography on small subsystems.Furthermore, we note that the size of the subsystems required for studying the entanglement scaling behavior does not depend on the overall lattice size.Using our simultaneous tomographic readout capability we are able to reconstruct the density matrix for various subsystems throughout our lattice by measuring the same number of Pauli strings necessary to reconstruct a single density matrix describing the largest desired subsystem.This feature makes superconducting qubit arrays desirable for studying entanglement within many-body systems.
Lastly, we provide an estimation of the measurement time needed to perform full tomography of a V = 8 subsystem.The measurement time is constrained by the time needed to upload measurement sequences from the measurement computer to the measurement hardware, plus the time needed to download measurement outcomes back to the measurement computer.In light of the rapid commercial development of improved measurement hardware, we expect these constraints to be assuaged in the near future.Based on the scaling shown in Fig. S18, a V = 8 subsystem requires 3×10 5 measurements for each of the 3 8 Pauli strings, or 2×10 9 total measurements.At the 100 µs measurement period (10 4 samples per second) used in this work, this corresponds to 55 hours of measurement time.Latency associated with measurement sequence upload can be reduced to constant time by programming the field programmable gate arrays onboard measurement hardware to synthesize combinations of pre-compiled pulse waveforms corresponding to each Pauli operator.Assuming state discrimination is performed onboard the measurement hardware, each measurement requires the transfer of 8 bits of data.Data transfer times will then be negligible, assuming 1 Gbit/s communication should become available.Adding extra time for occasional recalibration routines (to correct, for example, for slow drift in ⃗ Φ offset ), we, therefore, expect that near-term hardware developments will enable full tomography of a V = 8 subsystem in approximately 3 days.The Schmidt decomposition is a useful tool for representing entangled states.For a system AB in a pure state |ψ⟩, the theorem states that there exist orthonormal states |k A ⟩ for subsystem A, and |k B ⟩ for subsystem B, referred to as the Schmidt bases, such that where λ k coefficients are non-negative real numbers satisfying k λ 2 k = 1, and are known as Schmidt coefficients [23].The Schmidt coefficients provide insight into the "amount" of entanglement between the two subsystems.For two maximally entangled subsystems, each with n qubits, λ k = 1/ √ 2 n , which results in an entanglement entropy S 2 (ρ A ) = n.The pairwise entanglement for any two qubits in such a system is, therefore, exponentially small in the system size.In slightly entangled states, only a relatively small number of the Schmidt coefficients have a significant contribution to the weight of the state.Therefore, area-law states can be compressed and represented by a truncated Schmidt decomposition [22] |ψ trunc ⟩ =

S16. COHERENT-LIKE STATES IN 1D
In two dimensions, the hard-core Bose-Hubbard model is non-integrable, while in one dimension the model is integrable.While a general statement is beyond the scope of the present work, eigenstate thermalization is known to fail in many integrable systems.It is therefore interesting to compare the behavior of the two-dimensional (2D) system presented in this work with that of an analogous one-dimensional (1D) system.In this section, we present numerical simulations of coherent-like states prepared in a 14-qubit one-dimensional hard-core Bose-Hubbard chain, both with and without next-nearest neighbor exchange interactions.These simulations do not include decoherence.Unlike in two dimensions, in the one-dimensional system, the population does not reach a steady state in time of order 1/J.The average number of excitation continues to oscillate at for tens of characteristic periods 1/J, with the specifics of the oscillation amplitude and ringdown depending on δ (see Fig. S25a and S25b).The analogous dynamics in the presence of exaggerated next-nearest neighbor coupling J NNN = J/6 are shown in Fig. S26a and S26b.In Fig. S25d we present the two-point correlations averaged across all distance-3 qubit pairs in the 1D lattice as a function of time.Unlike in the 2D case, here the correlations fluctuate significantly in time, even after time 10/J.The lattice particle number ⟨n⟩ measured in experiments for states at approximately t ≈ 10/J.The number of excitations at equilibrium will vary depending on δ.When δ = 0, the lattice reaches half-filling, and as the magnitude of the detuning increases, the value of ⟨n⟩ at equilibrium decreases.

FIG. 1 .
FIG. 1. Experimental concept (a) Schematic for an example subsystem X of 4 qubits within a 16-qubit lattice.The subsystem has a volume of 4 (maroon sites) and an area of 8 (orange lines).(b) 2D hard-core Bose-Hubbard lattice emulated by the superconducting quantum circuit.Each site can be occupied by, at most, a single particle.(c) Energy E spectrum of the hard-core Bose-Hubbard lattice emulated by our device, shown in the rotating frame resonant with the lattice sites.The energy spectrum is partitioned into distinct subspaces defined by the total particle number n.(d) Scaling of the entanglement entropy S with subsystem volume V for an eigenstate at the center of the energy spectrum (orange line corresponding to the energy eigenstate highlighted by an orange oval in c) and an eigenstate at the edge of the energy spectrum (teal line corresponding to the energy eigenstate highlighted by a teal oval in c).(e) Change in the entanglement behavior, quantified by the geometric entropy ratio sV /sA, for states with n = 8. (f ) Schematic for the flip-chip sample consisting of 16 superconducting qubits.Optical images of the qubit tier (g) and the interposer tier (h) are illustrated with the qubits and the different signal lines false-colored.

FIG. 2 .
FIG. 2. Coherent-like state preparation.(a) Total number of particles ⟨n⟩ in the uniform lattice while driving the system on resonance for time t.Simulations do not include decoherence.The lattice reaches half-filling at equilibrium.The experiments are executed using the pulse sequence shown in the inset.(b) Probability of measuring different numbers of excitations in the lattice at three different times.The blue stars are from a Poisson fit to the excitation number distribution for Ω = J/2, with the dashed lines as guides to the eye.Simulated overlap of the prepared coherent-like state in steady state (t = 10/J) with drive strength Ω = J/2 and drive detuning δ = 0 (c), δ = 1J (d), and δ = 2J (e) with the hard-core Bose-Hubbard energy eigenstates.The different shades of red indicate the magnitude of the overlap between the prepared superposition states and energy eigenstates.Note the spectra are shown in the rotating frame of the lattice sites, not of the drive.(f ) Average 2-point correlator squared along the x-basis, |C x i,j | 2 , between qubit pairs at distance M for drive duration t = 10/J, strength Ω = J/2 and detuning δ from the lattice frequency.(g) Correlation length ξ x extracted using the 2-point correlators at different values of δ.

FIG. 3 .
FIG. 3. Entanglement scaling behavior (a) Examples of subsystems measured in our lattice.(b) The average subsystem entanglement entropy S2(ρX ) for subsystem with volume V for a state prepared with drive time t = 10/J, strength Ω = J/2 and detuning δ from the lattice frequency.(c) Subsystem entanglement entropy S2(ρX ) along dashed lines in panel (b) for detunings δ/J =0 (green), 0.9 (orange), and 2.1 (blue).Solid circles are experimental data averaged across all subsystems of each volume, with data for individual subsystems indicated by smaller, transparent circles.For δ/J =0 (green), open diamonds are Monte Carlo simulations of 20,000 samples indicating that deviations at subsystem sizes 5 and 6 arise partially from insufficient sampling.The gray region indicates the estimated classical entropy from dephasing (see Section S14 of the Supplementary Material for details).(d) The volume entanglement entropy, sV , and area entanglement entropy, sA, per site extracted using 163 different subsystems of various volumes and areas for states prepared with drive detuning δ.The error bars indicate ±1 standard error of the fit parameter.(e) The geometric entropy ratio sV /sA is used for quantifying the behavior of entanglement.States prepared with δ = 0 exhibit a strong volume-law scaling, and the states prepared with larger drive detuning values show a weaker volume-law scaling with an increasing area-law scaling.Dashed lines in (c-e) indicate results from numerical simulation.

FIG. 4 .
FIG. 4. Schmidt coefficients scaling across the spectrum.(a) The Schmidt coefficients for the decomposition of a five-qubit subsystem (highlighted with maroon color in the inset) and the remaining lattice.Coherent-like states are prepared with different drive detunings δ/J = 0, 0.9, 2.1.The coefficients in the left panel are calculated from experimental data and agree well with the simulated coefficients in the right panel.(b) The ratio of the largest and the k th -largest Schmidt coefficient λ 2 1 /λ 2 k , for k = 5, 10, 14, of coherent-like states prepared with different drive detunings δ in experiment.(c) The number of Schmidt coefficients required to represent the bipartition of the lattice with subsystem volume V = 3, 4, 5 with accuracy 1 − ϵ = 0.999.Each data point is averaged over all the measured subsystems of the same size, with the error bars indicating ±1 standard deviation.

S15
partment of Defense, or the Under Secretary of Defense for Research and Engineering.AUTHOR CONTRIBUTION AHK performed the experiments and analyzed the data with ITR.AHK and ITR performed the numerical simulations and analysis used for characterizing the system and validating the experimental results with assistance from ADP and YY.AHK and CNB developed the experiment control tools used in this work.AHK, CNB, ADP, LD, PMH, MH, and JAG contributed to the experimental setup.SEM, PMH, and MS designed the 4 × 4 qubit array.RD, DKK, and BMN fabricated the 4 × 4 flip-chip sample with coordination from KS, MES, and JLY.SG, JAG, and WDO provided experimental oversight and support.All authors contributed to the discussions of the results and to the development of the manuscript.S4 S3.Control and readout calibration S5 A. Flux crosstalk calibration S5 B. Time-domain flux pulse shaping and transient calibration S8 C. Pulse alignment S9 D. Readout optimization S10 E. State preparation and measurement (SPAM) error mitigation S12 F. Idling frequency layout S12 S4.Hamiltonian Characterization S13 A. Particle exchange strengths S13 B. Drive coupling amplitude S13 C. Drive coupling phase S14 S5.Hard-core Bose-Hubbard model energy spectrum S16 S6.Coherent-like states across the spectrum S17 S7.Impact of the drive coupling phase on experiments S18 S8.Measurement of global entanglement S19 S9.Tomography subsystems for studying the entanglement behavior S20 S10.Extracting s V and s A S22 S11.Measurement sampling statistics S23 S12.Extracting the entanglement scaling behavior in larger lattices S24 S13.Schmidt coefficients scaling S26 * karamlou@mit.edu† william.oliver@mit.eduarXiv:2306.02571v4[quant-ph] 25 Dec 2023 S2 S14.Estimating the impact of decoherence on measured system entropy S27 FIG. S1.Measurement setup.
FIG. S2.Learned Flux Crosstalk Matrices.The learned (a) DC flux crosstalk matrix and (b) fast flux crosstalk matrix for the 16-qubit device, each rescaled by a factor of 100 such that each element is a percentage.The fast flux microwave line for qubit 12 was broken, so we have no crosstalk information for qubit 12.
FIG. S3.Pulse alignment calibrations.(a)The timing for the Z pulse applied through each line is calibrated by sweeping the delay time of a Z pulse by tuning the qubit to the sweet spot while applying an XY pulse through the common resonator feedline at that frequency.(b) The timing for the XY pulse applied through each local line is similarly calibrated by applying a Z pulse with a constant delay characterized in (a) and sweeping the delay time of the XY pulse.
FIG. S4.Readout single-shot data.(a) We use an SVM to draw a discrimination line (dashed line) between the single-shot points prepared in the |0⟩ state (blue points) and |1⟩ state (red points) prior to readout.From this data, we obtain the readout fidelity FRO = 98.0%.(b) The same single-shot dataset is grouped into two clusters using a GMM with a cluster separation fidelity of Fsep = 99.9% FIG. S5.Readout optimization.(a) Experimentally measured readout optimization landscape using the readout fidelity.The plateau in the quantity is caused by qubit relaxation during readout.(b) Experimentally measured readout optimization landscape using the cluster separation fidelity.The star depicts the optimal parameters.(c) FRO and Fsep at different readout integration times.
FIG. S7.Histogram of coupling strengths among (a) nearest-neighboring sites and (b) next-nearest-neighboring sites.In a 4-by-4 square lattice, there are a total of 24 and 34 such couplings, respectively.
FIG. S8.Characterizing the drive coupling amplitudes.(a) Rabi oscillations observed by sweeping the drive pulse duration and amplitude.(b) Rabi rates extracted for different values of the pulse amplitude.(c) The amplitude of the relative drive coupling |αi| for different qubits in the lattice.
FIG. S9.Characterizing the drive coupling phase.(a) Pulse sequence used for interferometric pairwise relative drive coupling phase characterization.(b) Measured population on qubit "B" as a function of sweeping the exchange interaction time t and the relative drive phase ∆ϕ of the π 2 -pulse applied to qubit "B" before the readout.The gray dashed line indicates the relative common drive phase between the two qubits.(c) Population versus time for the 9 different 2 × 2 plaquettes; experimental data are shown as circles and the numerical simulations with the optimized drive coupling phases are shown as lines with colors corresponding to the qubit indices shown in parentheses.(d)The phase difference between the adjacent qubits extracted from the respective 2 × 2 plaquette fits.The two numbers in grey were unused because those bonds were captured through a higher-quality fit in the neighboring plaquette.The two red-shaded plaquettes were unused because of poor fit quality.
FIG.S10.Experimental data comparison to simulation.with (blue line) all phases are set to zero, (orange line) using phases characterized with the pairwise interferometry experiments, and (green line) with fine-tuned phase values.
FIG. S11.Energy spectrum comparison.The energy spectrum of the characterized system Hamiltonian (a), and the Hamiltonian without the next-nearest coupling (NNN) terms (b).The NNN couplings cause the spectrum to skew towards higher positive energies.(c) The difference between the magnitude of the highest eigenenergy Emax and the lowest eigenenergy Emin for the Hamiltonian with and without NNN terms.
FIG. S12.Entanglement scaling at different particle-number subspaces.(a) Subsystem entropy versus volume of states at the center and edge of the energy band for subspaces with n =5, 6, 7, and 8 particles.(b) Geometric entropy ratio sV /sA comparison for different particle-number subspaces.

x
FIG. S13.Impact of the drive coupling phase on correlations and entanglement.Simulations of (a) the average entanglement entropy of the 8-qubit subsystems, (b) the correlation lengths, and (c) the correlation matrices within the 4 × 4 lattice when driven on-resonance with and without the characterized drive phase.
FIG. S14.Entanglement formation time dynamics.(a),(b)The global entanglement build-up in the lattice after driving for time t with strength Ω = J/2 and detuning δ from the lattice frequency.Simulations do not include decoherence.

FIG
FIG. S15.Tomography subsystems.(a) Matching colors indicate the qubits in our lattice that have identical Pauli operators in each measured Pauli-string for tomography.(b) Distribution of the area (AX ) and the volume (VX ) of the measured subsystems.
FIG. S16.Visualizing entanglement in different subsystems.Scatter plots of the second Rényi entropy extracted of all 163 subsystems, shown as a function of the area (AX ) and the volume (VX ) of the each subsystem.Subsystem entropies are shown for the coherent-like states prepared at (left) δ = 0 and (right) δ/J = −2.1.
FIG. S18.Simulation of the measurement sampling statistics problem.The dark markers present the second Rényi entropy extracted from simulated tomography of subsystems as a function of the number of measurement samples of each Pauli string used for density matrix reconstruction.For each subsystem volume, results were averaged over the same subsystems used in the experiment.The probability distribution of measurement outcomes was determined from the simulated state prepared at Ω = J/2 and δ = 0.The light dashed lines represent the exact entropy of the simulated state at each volume.(a) Density matrices are reconstructed from the Stokes parameters obtained from Monte Carlo sampling of the probability distribution of eigenvalues of Pauli operator strings.(b) Density matrices are reconstructed using maximum-likelihood estimation on the bitstrings obtained from Monte Carlo sampling of the probability distribution of measurement outcomes of each Pauli operator string.Maximum-likelihood estimation is not used in (a), whereas the simulation in (b) uses the same reconstruction procedure that was used for experimental data.
FIG. S20.Scaling of the experimental approach in extracting the volume-to-area entropy per site ratio sV /sA.(a) The sV /sA ratio extracted for various states for 1d HCBH lattices with 16 and 25 qubits.(c) The sV /sA ratio variation for different superposition states as a function of the maximum measured subsystem size Vmax in 1d lattices.(c) The sV /sA ratio extracted for various states for 2d HCBH 4 × 4 and 5 × 5 lattices.(d) The sV /sA ratio variation for different superposition states as a function of the maximum measured subsystem size Vmax in 2d lattices.(e) The contrast of the sV /sA ratio extracted using maximum subsystem size Vmax between superposition states with r = 0.1 and r = 0.8.
FIG. S21.Schmidt decomposition (a)The 32 largest Schmidt coefficients for a bipartition of a volume-law eigenstate at the center (energy E ≈ 0) of the 16-qubit HCBH energy spectrum with an area-law eigenstate at the edge (energy E ≈ 15.2J) of the spectrum.

λ
FIG.S23.Entropy contribution from decoherence.An estimate of the entanglement entropy contribution to the measurement from spin-locking dephasing (with timescale 1/Γ).
FIG. S24.Time dynamics simulations.(a) Excitation population for driven 3 × 3 and 4 × 4 lattices.(b) The time evolution of the entropy of subsystems of different sizes in lattice driven with drive strength Ω = J/2 and detuning δ = 0. (c) The time dynamics of sV /sA ratios for various drive detunings δ.(d) The time dynamics of correlation lengths ξ x for various drive detunings δ.
FIG. S25.Simulations of coherent-like states in a one-dimensional 14 qubit system.(a) Measured total number of excitations ⟨n⟩ while driving the system with drive strength Ω = J/2 and various detunings δ for time t.(b) Total number of excitations versus δ, averaged between t/J = 15 and t/J = 20.(c) The second Rényi entropy of the subsystem consisting of the first seven qubits in the chain, shown as a function of time for various drive detunings.(d) The 2-point distance-3 correlator |C X M =3 | 2 , averaged over all distance-3 pairs in the lattice, as a function of time for various drive detunings.(e) The time-averaged 2-point correlator shown as a function of δ and pair distance M .The shaded region in a indicates the arbitrarily-chosen time window throughout which the results shown in subfigures b and e were averaged.
FIG. S26.Simulations of coherent-like states in a one-dimensional 14 qubit system with a next-nearest neighbor coupling of strength J NNN = J/6.(a) Measured total number of excitations ⟨n⟩ at drive strength Ω = J/2.(b) Total number of excitations versus δ, averaged between t/J = 15 and t/J = 20 (indicated by the shaded region in a).(c) The second Rényi entropy of the subsystem versus time, for the subsystem consisting of the first seven qubits in the chain.
FIG. S27.Extended data on generating coherent-like states.(a) Measured population on each qubit in the uniform lattice while driving the system on resonance for time t with different values of drive strength Ω. (b)The total number of excitation ⟨n⟩ in the uniform lattice while driving the system on resonance for time t.We observe excellent agreement between experimental data and numerical simulations using the characterized Hamiltonian parameters.

TABLE SI .
Active electrical components used in the experimental setup.

TABLE SII .
Sample parameters.We show the maximum transmon transition frequencies ω max q